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Explicit upper and lower bounds for the value F(u) of a linear functional F applied to a 
function v(x) defined on the interval 0<x<l are given when u is prescribed at the iV+1 
points i/N, i = 0, . . . , A r , and a bound for the integral of u [k]2 is known. These bounds 
are optimal in the sense that they are attained for functions satisfying the prescribed condi- 
tions. Their computation involves the inversion of a matrix of size k— 1 rather than N, 
which means that N is permitted to be large. 



1. Introduction 

Many problems in numerical analysis can be re- 
duced to approximating the value F(u) of a given 
linear functional F operating on an unknown element 
u of a linear vector space. The approximation is to 
be made in terms of a finite set of data concerning u. 
Thus, the values F x (u), . . . , F n {u) of N linear 
functionals acting on u may be given. For example, 
the F t (u) may be values of the function u at certain 
points Xi. If F(u) is the value of u at another point 
£, we have the problem of linear interpolation. If 
F(u) is an integral of u, we have the problem of nu- 
merical quadratures. If F(u) is the value of a deriv- 
ative, we have numerical differentiation. 

It was shown by M. Golomb and the author 2 that 
in order to obtain a finite interval in which the value 
F(u) must lie, one must be given the value of at least 
one nonlinear functional operating on u. The sim- 
plest case is that in which one is given a bound for a 
quadratic functional (u,u). In this case Golomb 
and the author 2 showed how to obtain the exact 
interval in which F(u) must lie when the values 
Fi(u), . . . , F N (u), and (u,u) are given. That is, 
upper and lower bounds for F(u) which are attained 
for some elements u satisfying the given conditions 
are found. The construction of these bounds re- 
quires the inversion of a matrix depending upon the 
functionals F, F h . . . , F N , and (u,u). 

In this paper we restrict our attention to a very 
simple case. We deal with a function u(x) of a 
single variable on the interval [0,1]. The given 
functionals F t are the values u(i/N) of u at the N-\-\ 
equally spaced points i/N, i=0, . . . , N. The 
quadratic functional is taken to be the integral of the 
square of the k th derivative of u. 

We think of the number of points JVas large, while 
the number k of the derivative will usually be small, 
say two or three. The matrix to be inverted is of 
size N. By making use of the equal spacing of our 



points, we shall reduce its inversion to that of a 
(k-l)X(k-l) matrix. Thus the problem of ob- 
taining best formulas for interpolation, quadratures, 
and numerical differentiation is made manageable 
even when the number of points involved becomes 
large. 

When F(u)= udx, our results yield as special 

Jo 

cases the best quadrature formulas of Sard 3 for 
k<3, N<6. 

2. Approximation Problem 

Let the values of the unknown function u(x) be 
given at the N+l evenly spaced points i/N, 
i=0, . . ., N. Let M 2 be a given bound for the 
integral of the kth derivative of u. 



(u,u) = 



J"- 

Jo 



™*dx<M 2 . 



(2.1) 



We assume that N>2k, k>2. 

Our problem is to approximate the value F(u) of a 
certain linear functional F applied to u. According 
to the theory in footnote 2, this is possible if and 
only if the functional F is bounded in the norm (2.1) 
for functions vanishing at the points i/N. That is, 
we must assume that there is a constant c such that 

F{v) 2 <c{ v [k]2 dx (2.2) 

for all k times differentiable functions v(x) such that 

Any linear combination of pointwise values or 
integrals of v and its derivatives up to order k-1 will 
satisfy this condition. 



1 Present address: University of Minnesota, Minneapolis, Minn. 

2 M . Golomb and II. F. Weinberger, Optimal approximations and error bounds, 
Symp. on Numer. Approx. (Univ. Wisconsin Press, Madison, Wis., 1959). 



3 A. Sard, Best approximate integration formulas; best approximation formulas. 
Am. J. Math. 71, 80-91 (1949). 
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To obtain a best estimate for F(u) we construct 
an auxiliary function u defined by the properties 



uW]( x )=0, 0<x<l, Nx?*0,l, . . ., N (2.4) 



and 



u"(0)=tt»'(l)=0 



i=0,l. 



.,N, 



l=k, . . ., 2k-2. (2.5) 



The function u and its first 2k-2 derivatives are 
continuous, while u {2k ~ l] is allowed to have jump 
discontinuities at the points i/N. 

The Green's function G(x; £) is defined by the 
properties 



0, 



d ik G 
da; 2 * 

<#0 



d'G 



dx 



T=0, 



0<X1, Nx?*N£, 0,1, ... N, 

i=0,l,. ,.,N, 

at x=0,l, l=k, . . . 2k— 2, 



d^fffc+O;*) d 2 *- 1 G f tt-0;£)_ 



dx 2k 



dx 2k 



-{-\f. 



y (2.6) 



Again G and its first 2k-2 derivatives are continuous, 
while the (2k — l)st derivative may have discon- 
tinuities at i/N. By integration by parts we find 
that 



«(«)=gflf<ft)«(^r). 



(2.7) 



where 

j/ i te)=(-ir i 

V-^^+Ojs) ^-'(7(^-0;^ 



da; 2 *- 1 



da; 2 ' 



(2i 



Also by integration by parts we have 

b k G(x; £) 



w(£)— w 



<«-r 



6 [fc]_^w 



dx A 



dx. (2.9) 



Applying the functional F and using Schwarz's 
inequality, we find that 

\f(u)-f(u)\ 2 < fv^-^ 1 " 1 ! 2 ^ r i^r^T^ir^ 

Jo J o I L ® x Jl 



-jj 



i.i*]- 



-uW\*dxFi{F,[G(l;;r,)]}. 



(2.10) 



(The symbol Frj [G (£; 77)] means that the functional 
F is applied to G considered as a function of rj for 



fixed £. F^{Fy,[G (£; 17)] } means that i^is then applied 
to the function F V [G]). We have used the property 



0ft; *)= 



r 






dz, (2.11) 



which follows from integration by parts. 
Another integration by parts shows that 



(2.12) 



P{#- u [fc] }dx=0. 

Hence, we can rewrite (2.10) as 

\F{u)-F{u)\ 2 <^M 2 -^u^ 2 dx^ 

(2.13) 

Once u and G are found, this inequality provides 
upper and lower bounds for F(v). These bounds are 
sharp in the sense that we can construct functions u 
satisfying (2.1) and having the given values u(i/N) 
for which the bounds for F(u) are attained. 
We write u in the form 



N+k-l 

i=-k+i 



(2.14) 



The centered difference operator A 2 * is defined by 

A2 ^=s(/+0 ( - i)t+lCi - (2 - i5) 

It is easily seen that A 2Jc \Nx— i\ 21i ~ l vanishes for 
\Nx—i\>k. Hence the sum in (2.14) has at most 
2k-l nonzero terms. 

Clearly the function (2.14) satisfies (2.4) and has 
the required continuity properties. 

The coefficients a t are to be determined by the 
conditons (2.5). Thus we must have 



N+k-l 

23 aA 2k \j-i\ 2k - 

i=-h+l 



-«(*> 



3=0,1, ...,N. 

(2.16) 



In order to apply the last line of (2.5), we first use 
partial summation to write 



5(a)=(— 1)* 2 D*a*D*|iVa5— *|*-', (2.17) 

i= — co 

where we have put a t =0 for i<—k, i>N+k, and 
where D k is the kth forward difference operator: 



D k c 



-s© 



(-lf +l c t+l 



(2.18) 



Since A*|iVx— i| a_1 =0 for [iVac — »| >ifc, (2.14) and 

(2.17) are independent of the values of a, for i < —k, 
i>N+k when 0<a;<l. We note that D k \Nx— ^l 2 * -1 
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is a polynomial of degree k—\ for Nx<i or Nx>i-\-k. 
Hence the List line of (2.5) will be satisfied if 



D k i 



0, 



-1, N-k+1, 



N—l. (2.19) 



The boundary value problem (2.5) is thus replaced 
by the system (2.16), (2.19) of N+2k—l linear 
equations in the N+2k— 1 unknowns a t . 

We turn to this problem of matrix inversion. 
Remark: We have assumed that k>2. The case 
k = \ is easily treated. Let u be the broken linear 
function coinciding with u at the points i/N and 
having its breaks at these points. Let G(x£) be 
the broken linear function with breaks at £ and the 
two neighboring points of the form i/N which 
vanishes at the points i/N, and whose derivative 
decreases by -1 at £. Then (2.13) gives optimal 
bounds for F(u) when k=\. 

3. Matrix Inversion 
We consider the system of linear equations (2.16), 



(2.19). Since A ? *|j- 



lepends only upon \j- 



and vanishes for \j~i\>k, (2.16) is a finite dif- 
ference equation or order 2k — 2. We solve it by 
means of a system of 2k — 2 independent solutions 
of the homogeneous equation. To find these solu- 
tions, we note that for any number z 



X) e iz A 2k \j-i\ 2k ~ 
f— M-l 



-2eV-U z (e z -l) 2k 



(£)" 



[(«•-!)- 



(3.1) 



forO<j<N. 

We define the polynomial in e z 



Q,(e g )m2(-l) l + l e-(e'-l) 



1+2 



\dz) 



[(e z -l)-'] 



(3.2) 



It is easily seen to be of degree I with leading 
coefficient 2. 

The Qi can be generated by the recursion 



Qo(y)=2 

The first few of these polynomials are 



Qi(y)=2(y+l), 

Q 2 (y)=2(tf+±y+l), 

Q 3 (y) = 2(tf+ntf+lly+l), 

Q i {y) = 2(y i +2Qy 3 -{- 66^+26?/+ 1) . J 



(3.3) 



(3.4) 



The coefficients of $2^-2 are the coefficients of the 
finite difference equation (2.16). It can be shown 



by induction that they are symmetric in the sense 
that 



(h (l)=y- l (My) 



(3.5) 



This means that the zeros of Q t occur in reciprocal 
pairs. When I is odd, one of the zeros is — 1. When 
I is even, Q can be written as a polynomial in (2/+?/ -1 ). 
Thus, the zeros of Q 2 k-2 can be found by solving an 
equation of degreed— 1 and a quadratic equation. 

It is shown by induction that the coefficients of 
Qi are positive and that the zeros of Q t are real and 
negative. The zeros of Q x separate thos of Qi+\. 

Let 2/i<C2/2<C • • • <d/2fc-2<C0 be the zeros of 

Q2 k- 2- 

#2*-2(2/„) = 0, V=l, . 



Then by (3.7) 



2/2*- l-v — 



2k-2. (3.6) 



(3.7) 



Because of (3.1) and (3.2), the functions a t =y v 
satisfy the homogeneous equations corresponding to 
(2.16). We now define 



*i= 




i<k-2. 
i>-k+2, 



(3.J 



It follows from the Lagrange interpolation 
formula [4] that the two definitions of \p t coincide for 
-k+2<i<k-2. Using (3.1) and (3.2), we find 
thai 



N+k-l 
i = -k-l 



^»|i-t|*-'=« 



Jp> 



j,p=0,... 9 N. (3.9) 



The same equation is satisfied by if/^p plus any 
linear combination of the y\. We add a linear 
combination such that the new function satisfies 
(2.19). Since \p t ^ p vanishes for i<p-\-k — 2 these con- 
ditions are already satisfied for i= — &+1, . . . , — 2. 
Therefore we add only functions which also satisfy 
these conditions. The Lagrange interpolation for- 
mula 4 shows that such linear combinations are 



furnished by the functions 7j i+a , a=l, 



Note that 
We let 



2k-2 yi-\-k-2 

" i= S (y,-i)*Q*-2(y r )' 



D k ni=4>i, i>-k+2. 



■*- ip Yi — p / 1 ^aVi 



, k, where 
(3.10) 

(3.11) 

(3.12) 



and determine the coefficients c a in such a way that 



* J. F. Stefiensen, Interpolation (Chelsea Press, New York, N. Y., 1950. 
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D k T ip =0 for i=-l, N-k + 1, . . ., 2V-1. This 
gives Cjfc=5 0p , and the equations 



k-i 



2 c a \// N _ k + a +p=D k \l/ N _ k+ p_ p 

a = l 

-topfN+e, P=l,...,k-1. (3.13) 

We let A a p be the inverse of the symmetric 

(k— l)X(k— 1) matrix \l/ N _ k+a+ p so that 5 



0=1 

(3.14) 
Then I\ p is given by 

/v-l 

a,0=l 

— i+1, • • ., N+lc-1, p=0, ...,N. (3.15) 
This function satisfies 

N+k-l 

1] r iP A*\i-j\*-i=8 jp , j,p=o,...,N 

(3.16) 
and 

Z>JIV=0, *=— ifc-hl, ..-,— 1, 

iV-fc+1, . . ., 2V-1; y=0, . . ., AT. (3.17) 
The solution of (2.16), (2.19) is given by 



Since k— 1 = 1, the matrix ^4 a/3 has the single element 
&J 1 . Moreover, 



Thus (3.17) gives 

V3(i-yr»)(i-yr^-^) y ;-^ 



(3.21) 



*=§r,,u(^} 



(3.18) 



Therefore r ip is the inverse matrix for the problem 

(2.16), (2.19). 

Our problem of matrix inversion has thus been 
reduced to the solution of polynomial equations of 
degree k—1 and two in order to find the y v and the 
inversion of the (k—1) -dimensional symmetric ma- 
trix l/^-fc+a+.tf. 

Example: We consider the case k=2. The zeros 
of Q 2 (y) are 

2/i=-2-V3, 

7/ 2 =-2+V3. (3.19) 

The function \[/ t defined by (3.10) is 

i<Q, 

*'-^ (*-<-,!) i>0. <320) 



» s The fact that tfw-jfc+a+$ is nonsingular follows from the uniqueness of u de- 
fined by (2.4), (2.5), and the linear independence of the functions A 2 k \Nx—i\* k ~ l 
and rn+a, 0=1, . . ., fe. 



r^H 



12(1— 7/r 2jV ) 



—i<i<y, l<y<iV— 1 



VS(l-y,-»)(l-y,- 



2(N-i)\ /) .p-i 



)yV 



12(l-y^) 

p<i<N+l, l<p<N-l 

(i-y[ 2i )y\- N 



i2(i-7/r^) 
i2(i-2/r 2iV ) 

4-V3-(4+V%f^ 



12(1-^) 



i<iV, 2?=iV 



i=0, p=0 



(i=— i, ^=o 

(i=iV+l, #=iv 

(3.22) 



4. Bounds 



We now return to the consideration of section 2. 
The function w is given by 



u(x)= S a^lNx-il*- 1 , 

|ATx— il</c 



(4.1) 



with 



>*=gr^(£). (4.2) 



Thus, we may write 

u(x) = ^gj(x)u(J^y 



where 



«— jfe+i 

|.ZV:r-i|<fc 



(4.3) 
(4.4) 



is the optimal approximation function corresponding 
to u(i/N) = 8 ij . 

In order to find the bounds (2.13) we need the 
integral of u [k]2 and 6r(x,£). Integrating by parts 
we find that 



[V , tfa=(«l)*gtt(^I5»-«W, (4.5) 



! 
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where [ ] j/N denotes the discontinuity in the function 
at j/N, with the convention 5 [2 *- 1] (0— )=u [2k ~ 1] 
(l + ) = 0. Applying partial summation to (4.1), 
we find that 

N+k-l 

u(x)= J2 INx-i^A^at. (4.6) 

i = -k+l 

Consequently, 

P uM 2 dx=(-l) k 2N 2k - 1 (2k-l)\ £ u (;^V%* 

= (-l) k 2N 2k - l (2k-l)\ 

gA^) u (i) A2n « (47) 

with the convention that 

A 2 /T„=Z?T^ for i=0, i=N. (4.S) 

We now construct the Green's function G(x£). 
We begin with a function H(x£) having the proper 
jump at x=£ and satisfying the condition that the 
derivatives of orders k, . . . 2k-2 vanish at the end 
points. Let p be any integer satisfying 



0<p<N£<p+2k-l<N. 



(4.9) 



It follows from the Lagrange interpolation formula 
(see footnote 4) that the function 



H(x;^)=\x^\^^ 

M =0 



N " 



"W«), (4.10) 



where 



&,(*)= 



( — l)"+ 1 JV**- 1 /2Jfc — 1 



-1)! V n ) 



(2k 

has the property 

H(*;f)s0 for z<-L a.>2±^zi. (4.12) 
Therefore, the function 



6,(351 * )= 2^t=0)l{ H(a!;a 



;=o 



0=0} 



(4.13) 

has all the properties (2.6). Thus, it is the Green's 
function. 6 We not that the sums in (4.4) and (4.13) 
involve at most 2k — 1 terms for each x and £. 

The bounds (2.13) for F(u) are now given ex- 
plicitly by 



N+k-l n / a \ 2 

F{u)- 2 Sr„F(A»|iVa!-t|»->)«(-^) 

i— *+i y=o VV/ 

1) ^ | /m*-2(-1)W 2 *- 1 (2A-1)! S A ? r ««(^)«(^) } 



< 



2(2fc- 



| F x [F ( (H(x; ?))l-.^5 + 1 i g r F(A*|ATa!-i|*-')F [ff (^J {)] } (4.14) 



where we again use the convention (4.8). If F is 
local in the sense that F(u) only involves the values 
of u in the neighborhood of a point, the sums in the 
second term on the right involve at most 2k— 1 values 
ofiandj. Example: Letk=2. r zj is given by (3.22). 
We find that 



A A \Nx—i\ : 



'2{4-6|A^-;| 2 +3|iVz-;| 3 } 
2{2-|A^x-t| 3 } 




\Nx—i\<l 



l<\Nx-i\<2 
\Nx-i\>2. 



(4.15) 



Thus the interpolation function #,(#) is given by 

g i (x) = 12(q+l-Nx)(Nx-q)[(q+2-Nx)T tti 
+ (Nx+l-q)T Q+h J ] + ( a +i-.Nx)% j ) 

+ (Nx-qY8 q+h j, 

where the integer q is defined by 

q<Nx<q+l. (4.17) 

If 1<N£<7V— 1, we let p be the integer such that 

p + l<N£<p + 2. (4.18) 



(4.16) 



6 Since Green's function is uniquely denned by (2.6), it does not depend upon 
the integer p used in the construction of II. The function Hdoes depend upon p. 
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Then we have 



3N d H(x; f ) 



(Ni- P )(N^-p-i)(p+2-m(p+s-m~ 





(Nx—pY 

m-p 

(Nx-p) 3 3(Nx-p-iy 
N^—p " N^—p—l 

(p+3-Nx) z 3(p+2-NxY 



p+s-m 

(p+3-NxY 

L p+z-W 



p+2-N% 



Nx<p } Nx>p+3 
p<Nx<p+l 

p+l<Nx<N£ 

N£<Nx<p+2 

p+2<Nx<p+3. 



(4.19) 



If 0<iVf< 1, the definition (4.18) gives p= — l, so 
that the function H defined by (4.19) does not vanish 
at x=0. In this case we simply subtract the 
function 



1 



(A^+l) _1 A 4 |iVx-l| 



(4.20) 



from the right-hand side of (4.19) with p= — l to 
obtain H(x,£). Similarly, if N—KNt<N, we 
subtract 

I (N+l-N^-'A^Nx-iN+l)^ (4.21) 
z 



from the right-hand side of (4.19) with p = N—2. 

Thus we can evaluate the bounds 4.14 explicitly. 
In the special case of linear interpolation we have 



If 



we find 



F(u)=u(fi. 



l<q<N?<q+l<N-l, 



(4.22) 



(4.23) 



TO(ff(z;|))]- g g T ti F(A'\Nx-i\>)F[H(jp {)] 



=N-\Nt-q)Kq+l-m) 1 -{ 4 



{' 



3(l-yr 2N ) 



[( q+ 2-NW(i-yr 2Q )(i-yi- 2(N - Q) ) 



+2(a+2-m)(N^-a+i)(i-yr 2a )(i-yr 2iN '^ 1) )yr 1 

+ (m-q+mi-yr 2{ « +1) )(l-y 1 - 2iN -«- 1) ) 



where 



2/i=— 2— V3. 



The first factor on the right of (4.14) does not 
approach zero as iV— > °° unless M 2 happens to be the 
exact value of the integral of u [k]2 . Thus the differ- 
ence between the best upper and lower bounds in 
the linear interpolation problem with iV+1 equally 

spaced points and with a given bound for I u" 2 dx 

is of the order N 3/2 . If a uniform bound for \u"\ 
is given, one can obtain bounds for u(£) which 
differ by a term of order TV -2 . This shows that a 



]}■ 



(4.24) 



bound for the square integral gives considerably 
less information than a bound for the maximum. 
The problem of finding best bounds when the maxi- 
mum of \u [k] \ is bounded is much more difficult than 
the problem treated here. 
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paper. 
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